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ABSTRACT 


The  two-dimensional,  steady,  compressible,  inviscid  flow  of 
ionized  gas  through  linear,  segmented  electrode,  magnetohydrodynamic 
(MHD)  channels  is  computed  in  the  plane  of  the  applied  electric  field. 
The  solutions  obtained  for  the  gas  dynamic  and  electrical  quantities 
satisfy  the  coupled  fluid  mechanical  conservation  laws  and  Maxwell's 
equations  at  all  points  of  the  flow  field  simultaneously.  The  nonuniform 
profiles  which  are  obtained  include  the  effects  of  transverse  variations 
in  incoming  gas  dynamic  profiles,  local  Hall  currents  in  Faraday  chan¬ 
nels,  transverse  currents  in  Hall  devices,  and  axial  magnetic  field 
gradients . 
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SECTION  I 
INTRODUCTION 


The  channel  flow  of  an  electrically  conducting  gas  has  been  of 
engineering  interest  in  recent  years  because  of  the  capability  of  pro¬ 
ducing  an  electromagnetic  interaction  with  the  fluid  which  results  in 
either  energy  addition  to  or  extraction  from  the  flow,  thus  causing  an 
acceleration  of  the  fluid  or  a  generation  of  electrical  power.  In  order 
to  make  theoretical  predictions  of  the  performance  of  these  devices, 
previous  investigators  introduced  assumptions  to  obtain  tractable  formu¬ 
lations  which  led  to  analytical  models  at  some  variance  to  their  physical 
counterparts.  Representative  of  these  approaches  are  the  solutions 
discussed  in  Ref.  1.  These  include  exact  solutions  to  certain  simpli¬ 
fied,  linearized  versions  of  the  magnetohydrodynamic  (MHD)  equations 
(Chapter  10  of  Ref.  1),  and  a  discussion  of  numerical  solutions  to 
approximations  of  the  describing  equations  (Chapter  11  of  Ref.  1).  The 
numerical  cases  considered  there  utilize  the  quasi -one -dimensional 
approximation  to  MHD  channel  flow.  Another, type  of  analytical  model 
which  has  been  studied  is  exemplified  by  Ref.  2  in  which  the  two- 
dimensional,  compressible,  laminar  boundary  layer  form  of  the  equa¬ 
tions  of  motion  are  solved  across  the  entire  channel  in  the  E -plane. 
There  the  assumptions  are  made  that  the  gas  obeys  an  ideal  equation  of 
state,  the  insulator  walls  are  parallel,  the  ratio  of  the  electric  field 
components  Ey/Ex  is  constant,  and  Ex  is  a  function  only  of  the  x  co¬ 
ordinate.  These  latter  two  assumptions  remove  the  requirement  for 
solving  an  elliptic  equation  for  the  electrical  properties.  The  analysis 
has  been  formulated  in  a  very  general  manner  in  Ref.  3  where  solutions 
have  been  obtained  that  include  finite  chemical  reaction  rates,  elec¬ 
trical  thermal  nonequilibrium,  thermal  and  concentration  diffusion, 
electron  energy  relaxation,  turbulent  boundary  layers,  and  a  one¬ 
dimensional  inviscid  core.  The  electrical  portion  of  the  problem  is 
solved  by  combining  separate  solutions  for  an  electric  current  stream 
function  in  the  inlet  region,  exit  region,  and  periodic  main  part  of 
Faraday  accelerator  channels.  The  geometry  considered  permits 
diverged  electrode  walls,  and  assumes  parallel  insulator  walls. 

The  analysis  performed  in  this  investigation  considers  compress¬ 
ible,  inviscid,  two-dimensional  flow  between  the  electrode  walls.  The 
distinguishing  feature  of  this  study  is  that  the  calculation  of  the  coupled 
gasdynamic  and  electrical  variables  is  performed  in  a  manner  such 
that  the  fluid  mechanical  conservation  laws  and  Maxwell's  electrical 
differential  equations  are  satisfied  simultaneously  throughout  the  entire 
field.  In  particular,  this  means  that  the  solenoidal  property  of  current 
density  and  the  irrotational  property  of  the  electric  field  are  satisfied 
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at  all  points  of  the  flow  field  where  electric  current  is  flowing.  The 
significance  of  this  is  that  it  is  possible  to  observe  the  changing  char¬ 
acter  of  the  current  paths  down  the  channel,  which  would  not  be  possible 
if  separate  solutions  to  the  electrical  portion  of  the  problem  were  ob¬ 
tained  over  portions  of  the  flow  field  and  then  joined  together.  The 
results  of  this  study  show  that  there  exists  a  nonnegligible  coupling  in 
the  currents  paths  along  the  channel.  The  mass  and  current  continuity 
equations  are  formulated  in  a  manner  which  permits  the  B-walls,  as 
well  as  the  E -walls,  to  be  nonparallel. 

In  order  to  demonstrate  the  physically  realistic  effects  of  the  MHD 
variables  being  coupled  together  at  all  points  of  the  flow  field,  the 
elliptic  differential  equations  for  electrical  variables  are  solved 
throughout  the  portions  of  the  channels  where  the  magnetic  field  is 
nonzero  (the  B- field  is  assumed  to  extend  beyond  the  electrically 
powered  portions  of  the  channels).  Attempts  to  achieve  the  numerical 
solutions  of  these  elliptic  equations  by  iterative  methods  were  unsuc¬ 
cessful.  This  was  primarily  a  consequence  of  the  large  gradients  in  the 
unknown  variables  and,  secondarily,  a  result  of  mixed  boundary  condi¬ 
tions  in  the  Hall  channel  analysis.  It  was  only  when  a  direct  nonitera¬ 
tive  method  given  in  Ref.  4  was  used  that  the  solutions  were  effected. 

The  results  of  the  numerical  experiments  performed  in  the  course  of 
this  study  indicate  that  this  method  can  be  rendered  stable  on  a  grid- 
work  limited  in  size  only  by  available  computer  memory,  in  contrast 
to  other  direct  methods  which  are  subject  to  truncation  error  instabili¬ 
ties  when  the  system  is  sufficiently  large.  This  particular  method 
seems  to  provide  not  only  a  sine  qua  non  for  problems  with  numerical 
complexities  as  severe  as  those  of  this  study,  but  a  more  efficient  tech¬ 
nique  for  simpler  problems  than  iterative  methods  provide. 

The  work  described  in  this  report  was  carried  out  with  the  follow¬ 
ing  objectives:  (1)  to  develop  computation  techniques  for  two-dimensional 
MHD  channel  flow,  (2)  to  assess  the  utility  of  using  quasi-one-dimensional 
calculations  for  design  purposes,  and  (3)  to  test  the  capability  of  the 
theory  to  predict  results  which  will  compare  favorably  with  available 
experimental  data.  The  second  objective  is  tested  by  averaging  the  two- 
dimensional  results  in  a  direction  transverse  to  the  channel  centerline, 
and  then  comparing  these  values  with  those  obtained  from  corresponding 
uniform  profile  solutions.  The  correspondence  of  the  two  methods  of 
solution  is  established  by  matching  the  average  values  of  the  assumed 
two-dimensional  profiles  at  the  channel  entrance  to  the  specified  one¬ 
dimensional  variables  characterizing  the  flow. 

The  third  objective  is  tested  by  comparing  theoretical  predictions 
with  available  impact  pressure  measurements  made  at  the  exit  of  an 
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accelerator  channel.  Since  the  initial  profiles  entering  the  channel  were 
not  measured,  an  empirical  method  is  employed  in  choosing  the  profiles 
that  gives  agreement  with  the  experimental  impact  pressure  measured 
in  the  B-field  plane  for  a  given  accelerator  run.  The  details  of  the 
method  are  described  in  Ref.  5. 


SECTION  II 

DESCRIBING  EQUATIONS  AND  ASSUMPTIONS 


The  analytical  model  of  the  MHD  flow  field  considered  is  obtained 
by  introducing  simplifying  assumptions  into  the  equations  that  express 
the  conservation  of  mass,  momentum,  and  energy  of  an  ionized  gas 
flowing  through  crossed  electric  and  magnetic  fields.  The  solutions  to 
these  equations  describe  the  flow  in  the  plane  of  the  applied  electric 
field  (the  x-y  plane),  while  at  the  same  time  they  include  the  effect  of 
diverging  conductor  walls  since,  in  general,  the  channel  dimension  W 
in  the  z  direction  is  both  finite  and  a  function  of  x.  Thus,  both  dimen¬ 
sions,  height  and  width,  are  permitted  to  vary  with  x.  In  order  to 
analyze  such  a  three-dimensional  physical  model  with  two-dimensional 
equations,  the  assumption  of  source  flow  is  made  (this  is  described  in 
detail  in  Appendix  II).  This  is  supplemented  with  the  additional  gas 
dynamic  assumptions  that  the  flow  is  steady,  inviscid,  and  nonheat¬ 
conducting.  It  is  further  postulated  that  the  magnetic  field  is  applied 
in  the  z  direction  and  that  the  induced  magnetic  field  is  negligible;  the 
electric  field  and  current  density  vectors  have  components  only  in  the 
x  and  y  directions.  The  formulas  for  the  variables  which  enter  into 
the  describing  equations  are  discussed  in  the  following  sections. 


2.1  GAS  DYNAMIC  FORMULAS 


The  fluid  mechanical  relationships  are  derived  from  the  conserva¬ 
tion  of: 

Mass  (Source  Flow) 


Momentum 


Energy 


7  .  WpV  =  0 


DV  -  - 

p  —  +  vp  =  jxB 
Dt 


P 


DH 

Dt 


(1) 

(2) 

(3) 


3 


AEOC-TR-71-181 


The  analysis  is  performed  using  a  streamtube  method.  For  this 
purpose  a  transformation  of  coordinates  is  introduced  which  uses  the 
von  Mises  variables  (x,  ip)  where  is  a  normalized  stream  function 
defined  to  satisfy  continuity  of  mass  and  possessing  the  properties 
that  p  is  constant  along  streamlines,  and  the  increment  in  ip  between 
streamlines  is  directly  proportional  to  the  mass  flow  in  the  streamtubes 
they  bound.  The  system  of  gas  dynamic  equations  is  closed  by  postu¬ 
lating  that  the  streamtube  slopes  in  the  x-y  plane  vary  linearly  in  the 
transverse  direction  from  the  centerline  to  the  walls 

dy(x,i|f)  =  y  d  £ 

dx  d  dx  2  (4) 

2 


Although  the  problem  is  solved  using  von  Mises  coordinates, 
several  of  the  describing  equations  can  be  written  in  a  more  compact 
form  by  using  natural  coordinates  in  the  notation.  Thus,  the  use  of  s 
and  n  as  subscripts  indicates  vector  components  in  these  directions, 
and  gradients  in  these  directions  can  be  related  to  gradients  in  the 
von  Mises  coordinates  from  geometrical  considerations. 


When  the  assumption  of  steady  flow  is  incorporated  into  Eqs.  (1) 
to  (3)  the  conservation  relations  are 


Mass  (Source  Flow) 

d*  2W 

—  =  —  pq 

dn  m 

Momentum 

d  q2  dp  dx 

p  - - —  +  t—  =  (-r— )  j  B 

dx  2  dx  ds  n 

dv  dt)  dp 

PU  d^  +  d7  d7  =  ’  dxB 

Energy 

^  .  q2  -*  -* 

pu  —  (h  +  =  E*j 


(5) 


(6) 

(7) 


(8) 


Subsequent  to  the  solution  of  these  equations  the  streamline  locations 
are  determined  from  the  integrated  form  of  the  definition  of  the  mass 
stream  function 

/t 

dij;  d 

pu  ”  2  (9) 

-1 
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In  order  to  facilitate  the  numerical  integration  of  the  equations  of 
motion,  they  are  transformed  into  an  equivalent  system  of  equations  as 
explained  in  Appendix  III. 

Separate  calculations  were  performed  where  the  inertial  term 

0V 

pu  in  the  transverse  momentum  relationship,  Eq.  (7),  initially  was 

included  in  the  pressure  prediction  scheme  and  subsequently  was  not  in¬ 
cluded.  The  results  of  the  calculations  were  virtually  insensitive  to  the 
choice  of  these  two  different  assumptions,  and  hence  all  the  results  of 
this  report  were  produced  neglecting  this  term  since  it  simplifies  the 
numerical  analysis. 

To  complete  the  formulation  of  the  gas  dynamic  flow,  relations  are 
needed  to  compute  the  thermodynamic  and  the  electrical  transport 
properties  of  the  working  fluid.  The  accelerators  considered  in  this 
report  use  either  air  or  nitrogen.  The  generator  is  supplied  with  the 
products  of  combustion  gases  (toluene  plus  methanol;  7C7Hg  +  2  CH3  OH  + 
66  O2  +  67  N2). 

The  thermodynamic  functions  derived  for  air  and  nitrogen  (N2) 
assume  these  gases  to  be  in  chemical  equilibrium  and  include  real-gas 
effects.  They  have  been  derived  by  postulating  functional  relationships 
which  simultaneously  satisfy  a  thermal  equation  of  state  and  the  second 
law  of  thermodynamics.  Additional  considerations  involved  in  the  selec¬ 
tion  of  these  relationships  are  that  they  should  provide  accurate  curve 
fits  to  the  exact  values  given  in  Refs.  6  and  7.  The  detailed  results  are 
presented  in  Appendix  IV. 


The  thermodynamic  functions  of  combustion  gas  products  are  ob¬ 
tained  from  curve  fits  to  unpublished  data  computed  at  AEDC. 


The  electrical  conductivity  of  seeded  air  and  nitrogen  is  computed 
from  the  method  described  in  Ref.  8.  The  method  incorporates  the 
known  theoretical  expressions  for  electrical  conductivity  of  a  singly 
ionized  plasma  in  the  two  limiting  cases  of  fully  ionized  and  slightly 
ionized  plasma  into  a  form  proposed  by  Lin,  Resler,  and  Kantrowitz 
(Ref.  9)  to  approximate  the  conductivity  (denoted  by  aQ  in  this  approxi¬ 
mation)  between  these  two  limiting  cases 


+  o 


(10) 


where  C7en  is  determined  by  electron-neutral  collisions  and  cjej  by 
electron-ion  collisions.  Subsequently,  Demetriades  and  Argyropolous 
(Ref.  10)  have  developed  a  higher  order  approximation,  written  as 
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CJ 


(11) 


in  the  notation  of  Ref.  8.  Equation  (11)  is  found  by  using  Grad's 
13 -moment  approximation  to  derive  a  generalized  Ohm's  law  for  a 
multicomponent,  nonisothermal  plasma,  including  temperature  and 
pressure  gradients,  as  a  function  of  the  electric  and  magnetic  fields. 
The  electrical  conductivity  is  obtained  in  two  successive  approximations 
which  are  equivalent  to  the  first  and  second  approximations  in  the  ex¬ 
pansion  of  Sonine  polynomials  in  the  Chapman-Enskog  approach. 


The  Hall  parameter  is  computed  from 


U)T 


(12) 


where  ene  is  the  product  of  electronic  charge  times  electron  number 
density. 


The  electrical  conductivity  and  Hall  parameter  of  a  seeded  com 
bustion  gas  product  are  obtained  from  curve  fits  to  unpublished  data 
computed  at  AEDC. 


2.2  MAGNETIC  FIELD 

It  is  assumed  that  the  magnetic  Reynolds  number  is  sufficiently  low 
so  that  no  magnetic  field  is  induced  by  the  flow  of  electrical  current. 

It  is  further  assumed  that  the  only  nonvanishing  component  of  the 
applied  magnetic  field  vector  is  in  the  z  direction,  and  that  this  com¬ 
ponent  is  a  function  of  the  x  coordinate  only.  A  functional  relationship 
B(x)  is  determined  for  each  channel  by  approximating  the  measured 
distribution  with  an  analytical  expression. 


2.3  ELECTRIC  FIELD  AND  ELECTRIC  CURRENT  DENSITY 

The  method  for  posing  the  analytical  model  describing  the  elec¬ 
trical  variables  in  a  given  channel  is  determined  by  the  manner  in 
which  the  boundary  conditions  are  prescribed.  In  the  case  of  the 
Faraday  channels  the  normal  component  of  current  density  is  specified 
on  all  the  boundaries hence,  it  is  natural  to  solve  a  partial  differential 
equation  for  a  current  stream  function  since  in  this  formulation  the 
boundary  conditions  possess  the  desideratum  of  being  of  the  Dirichlet 


6 


AEDC-TR-71-181 


type.  The  formulation  of  the  current  stream  function  equation  is  per¬ 
formed  in  Appendix  V  and  the  method  of  solution  is  described  in 
Appendix  VII.  In  the  case  of  the  Hall  channel,  the  terminals  are  equi- 
potential  surfaces;  it  is  natural  to  solve  a  differential  equation  for  an 
electric  potential  in  this  application.  The  formulation  of  such  an  equa¬ 
tion  is  performed  in  Appendix  VI,  and  the  method  of  solution  is  described 
in  Appendix  VII. 


Subsequent  to  the  solution  of  the  appropriate  elliptic  equation  for  ? 
or  ,  the  electric  field  and  current  density  are  related  through  an  Ohm's 
law.  The  formulas  which  are  combined  to  describe  the  electrical  por¬ 
tion  of  the  problem  are  Faraday's  law 

V  x  E  =  0  <13) 


continuity  of  current  in  two-dimensional  source  flow  (this  is  derived  in 
Appendix  II) 


and  an  Ohm's  law 


n 


7  •  Vfj  =  0 


a 

1  +  (idt)2 

q 

1  +  (out)2 


[Es  -  WT(En  -  qB)] 
[En  -  qB  +  U)T Eg  ] 


(14) 


(15) 


Equation  (15)  implies  that  the  total1  energy  addition  occurring  in  the 
energy  Eq.  (8)  can  be  written  as 

i  a  +j  2 

E-j  =  -  +  q  jn  B  (16) 

CT 


In  all  cases,  the  accelerator  cathodes  are  assumed  to  be  contained 
in  the  upper  walls. 


SECTION  III 
RESULTS 


A  series  of  flow-field  profiles  has  been  computed  in  order  to 
demonstrate  quantitatively  the  two-dimensional  characteristics  of  the 
variables  describing  MHD  channel  flow,  as  well  as  to  show  the  differ¬ 
ences  between  two-dimensional  and  quasi-one-dimensional  assumptions. 
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The  two-dimensional  characteristics  are  displayed  by  assuming  shapes 
for  the  initial  profiles  of  temperature  and  velocity,  and  then  following 
the  development  of  these  profiles  through  the  length  of  the  channel.  In 
addition,  the  integration  of  the  equations  of  motion  concurrently  with 
Maxwell's  equations  gives  the  current  components  in  the  x-y  plane  that 
are  induced  by  the  interaction  of  the  gas  dynamic  and  electromagnetic 
fields. 

The  theoretical  calculations  have  been  carried  out  for  flows  through 
three  accelerators  and  one  Hall  generator  that  are  analytical  models 
simulating  channels  constructed  at  AEDC.  These  will  be  referred  to  as 
Accelerator  A,  Accelerator  B  (a  20-MW  accelerator  designed  for  oper¬ 
ation  in  the  AEDC  LORHO  Pilot  Complex;  the  actual  slant-wall  boundary 
conditions  were  approximated  by  those  appropriate  to  a  Faraday  ac¬ 
celerator),  and  Accelerator  C  (described  in  Ref.  11).  The  generator 
calculations  were  performed  for  a  linear  (square  cross  section)  device, 
as  an  approximate  model  of  the  circular  cross  section  channel  described 
in  Ref.  12.  The  assumption  of  linearity  was  made  because  the  inter¬ 
action  of  the  applied  magnetic  field  with  the  flow  in  the  circular  chan¬ 
nel  produces  three-dimensional  effects. 


3.1  FARADAY  ACCELERATORS 

In  order  to  show  the  differences  in  accelerator  results  obtained  by 
solving  the  equations  of  motion  in  one  and  two  dimensions,  the  two- 
dimensional  profiles  are  averaged  in  the  transverse  direction  and  then 
one -dimensional  solutions  are  computed,  using  in  both  instances  the 
same  channel  entrance  values  of  appropriately  selected  variables  which 
characterize  the  flow.  In  particular,  the  following  initial  values  were 
matched: 


Accelerator  A 

Static  pressure 

Mass  averaged  total  enthalpy 

Mach  number 

Accelerator  B 

Static  pressure 

Mass  averaged  total  enthalpy 

Mass  flow 

The  two-dimensional  computations  incorporate  the  above  matching  con¬ 
ditions  in  conjunction  with  additional  assumed  conditions  at  the  channel 
entrance  as  given  in  Tables  I  and  II  (Appendix  VIII).  For  each  Accel¬ 
erator  A  and  B,  two  cases  of  two-dimensional  calculations  have  been 
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performed,  one  in  which  the  initial  temperature  and  velocity  profiles 
are  uniform  in  the  transverse  direction  and  a  second  in  which  the  initial 
profiles  are  nonuniform.  A  characteristic  feature  of  the  temperature 
profiles  is  that  initially  uniform  profiles  develop  overshoots  at  the  walls 
and  initially  nonuniform  profiles  tend  to  flatten;  the  latter  effect  is 
particularly  significant  and  is  primarily  attributable  to  the  joule  heating 
that  occurs  at  the  electrode -insulator  junctions.  These  results  are 
shown  in  Figs.  1  and  2  (Appendix  I).  The  differences  in  the  electric 
current  paths  are  shown  in  Figs,  lc  and  2c  (resulting  from  initially 
uniform  gas  dynamic  profiles),  and  Figs.  Id  and  2d  (resulting  from 
initially  nonuniform  profiles).  .  Two  significant  characteristics  of  the 
electric  current  paths  are  found  to  be  ( 1)  appreciable  eddy  current  flow 
both  upstream  and  downstream  of  the  pow'ered  section  -  in  fact,  the 
length  of  the  current  carrying  zone  is  on  the  order  of  twice  the  length 
of  the  powered  section,  and  (2)  the  current  path  between  the  ends  of  the 
powered  section  is  directed  approximately  straight  across  the  channel 
if  the  transverse  temperature  profile  is  nearly  uniform,  but  it  follows 
an  irregular  path  (jx  is  comparable  in  magnitude  to  jy)  if  the  tempera¬ 
ture  profile  is  nonuniform.  In  all  the  test  cases  that  were  computed  for 
these  two  accelerators,  it  was  found  that  the  wall  temperature  at  the 
downstream  electrode-insulator  junctions  is  a  highly  nonlinear  function 
of  the  electrode  current  density  at  the  points,  which  in  fact  determined 
the  downstream  values  of  jn  that  were  specified  in  the  analysis. 

The  typical  differences  in  downstream  values  between  variables 
calculated  from  (1)  quasi-one-dimensional  analyses  and  (2)  transverse 
averages  of  two-dimensional  analyses  initiated  with  nonuniform  pro¬ 
files  are  as  follows: 


p  9  percent 

u  2  percent 

T  7  percent 

H  3  percent 

Calculations  have  been  performed  for  flow  through  Accelerator  C 
to  test  the  capability  of  the  theory  to  predict  results  which  will  compare 
favorably  with  available  experimental  data.  This  objective  is  tested  by 
comparing  theoretical  predictions  with  impact  pressure  measurements 
made  at  the  channel  exit.  Since  the  initial  profiles  entering  the  channel 
were  not  measured,  an  empirical  method  is  employed  in  choosing  the 
profiles  that  gives  agreement  with  the  experimental  impact  pressure 
measured  in  the  B-field  plane  for  a  given  run.  The  details  of  the  method 
are  described  in  Ref.  5.  It  contains  the  premise  that  since  the  de¬ 
scribing  equations  are  inviscid  the  initial  profiles  should  include  non¬ 
uniformity  to  compensate  for  the  effects  of  boundary -layer  buildup. 
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The  salient  feature  of  the  technique  is  the  hypothesis  that  the  appropriate 

8u 

amount  of  initial  nonuniformity  is  that  which  results  in  -5—  =  0  in  the 

9z  wall 

accelerator  exit  region.  In  order  to  complete  the  specification  of  the 
initial-boundary  value  problem  it  is  necessary  to  prescribe  the  applied 
current  on  the  electrode  walls.  In  the  analyses  for  Accelerators  A 
and  B  it  is  stipulated  that  jn(x)  is  the  same  on  the  cathode  and  anode 
walls.  This  is  done  in  view  of  the  fact  that  the  segmented  electrodes  in 
those  channels  are  small  in  length  compared  to  channel  height.  The 
postulated  boundary  condition  that  jn  is  a  specified  function  of  x  would 
seem  appropriate  since  this  implies  infinitely  fine  segmentation.  This 
same  assumption  was  used  in  an  initial  attempt  to  calculate  the  flow 
through  Accelerator  C  even  though  this  device  has  large  electrodes.  The 
solutions  for  the  gas  dynamic  variables  were  almost  symmetric  about 
the  centerline  (Fig.  3)  in  contrast  to  experimental  results  (impact  pres¬ 
sures  measured  in  the  exit  region)  for  this  channel  (Fig.  3c),  which 
demonstrates  the  necessity  of  including  current  concentration.  Subse¬ 
quent  calculations  have  been  made  in  which  the  specified  jn  distribution 
on  the  walls  has  been  chosen  to  cause  the  current  to  flow  in  a  direction 
generally  from  the  anode  downstream  corner  to  the  cathode  upstream 
corner,  which  is  known  to  occur  from  both  theoretical  considerations 
and  physical  observations.  The  results  for  impact  pressures  (Fig.  3c)i 
do  indicate  asymmetry,  and  furthermore,  the  larger  theoretical  values 
of  impact  pressures  occur  along  the  anode  wall  in  accord  with  experi¬ 
ment.  However,  although  the  average  impact  pressures  show  favorable 
agreement  between  theory  and  experiment,  the  theoretical  results  do  not 
display  nearly  as  much  nonuniformity  as  the  measured  quantities.  The 
current  paths  in  the  powered  portion  of  the  accelerator  (Fig.  3d)  naturally 
reflect  the  assumed  boundary  constraints  on  (jn)wan*  but  •‘•n  t'le  entrance 
and  exit  regions  the  fringing  that  occurs  is  relatively  independent  of  the 
assumption  that  is  made  of  current  distribution  along  the  walls.  Thfcu'e 
theoretical  results  were  obtained  by  assuming  that  the  magnitude  of 
electric  current  applied  to  the  channel  was  that  which  had  resulted 
experimentally  when  the  applied  voltage  equaled  four  times  the  induced 
voltage.  Another  attempt  to  compare  theory  and  experiment  was  per¬ 
formed  at  a  higher  power  level  (applied  voltage  equal  to  six  times  in¬ 
duced  voltage)  but  resulted  in  numerical  instabilities  that  could  not  be 
resolved  as  a  consequence  of  the  temperature  becoming  unbounded  on 
the  anode  wall. 

3.2  HALL  GENERATOR 

The  one-  and  two-dimensional  Hall  generator  programs  have  been 
used  to  compute  a  case  simulating  a  design  point  of  the  LORHO  20-MW 
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generator  described  in  Ref.  12.  The  values  of  the  gas  dynamic  state 
upstream  of  the  magnetic  field,  terminal  voltage,  and  seed  rate  are  the 
same  for  both  methods.  It  is  found  that  the  two-dimensional  solution 
has  higher  values  of  pressure,  temperature,  and  axial  current  (Fig.  4). 
Since  the  potential  difference  between  the  terminals  is  the  same  for 
both  analyses  the  higher  current  implies  a  greater  generation  of  power 
(greater  by  16  percent  for  the  case  selected).  The  change  in  pressure 
across  the  channel  is  also  found  to  be  16  percent.  The  electrical  cur¬ 
rent  paths  and  Hall  voltage  distribution  are  shown  in  Figs.  4c  and  d, 
respectively.  The  current  paths  are  the  resultant  of  the  generated  axial 
current  and  the  transverse  current  which  is  short-circuited  through  the 
external  circuitry.  The  values  of  Hall  voltage  obtained  from  one-  and 
two-dimensional  approximations  to  the  describing  equations  compare 
closely  (Fig.  4d),  differing  by  3  percent  at  most. 


SECTION  IV 
CONCLUSIONS 


The  comparison  of  the  averaged  two-dimensional  with  the  quasi- 
one-dimensional  calculations  of  gas  dynamic  variables  shown  in 
Figs,  lb  and  2b  indicates  that  the  closeness  of  agreement  depends  upon 
whether  the  gas  influx  to  the  channel  in  question  is  uniform  (using  com¬ 
pression  heated  flow  as  a  gas  source,  for  example)  or  nonuniform  (such 
as  that  emanating  from  an  electric  arc  heater).  In  the  former  case  the 
comparison  is  quite  close,  which  would  seem  to  indicate  that  in  such  an 
application  the  quasi -one -dimensional  solution  could  be  used  for  the 
purpose  of  carrying  out  design  calculations.  In  the  latter  case  (initial 
nonuniform  profiles)  the  differences  can  be  as  much  as  10  percent, 
which  would  tend  to  indicate  a  greater  need  to  use  two-dimensional 
solutions  for  performance  predictions  in  these  applications.  In  both 
instances  there  is  a  dissipation  of  applied  electric  power  because  of 
current  fringing  in  the  end  regions,  the  effective  current  carrying  zone 
being  almost  twice  as  long  as  the  length  of  the  electrode  section. 
Initially,  nonuniform  profiles  result  in  additional  power  dissipation  in 
the  electrode  region  as  is  shown  by  the  existence  of  appreciable  local 
Hall  currents  (on  the  order  of  the  applied  Faraday  current)  even  though 
the  average  Hall  current  is  constrained  to  be  zero.  In  any  case  if  the 
flow  at  the  accelerator  exit  is  to  be  used  as  a  uniform  test  medium,  a 
two-dimensional  solution  will  have  to  be  carried  out  for  the  particular 
inlet  condition  considered.  The  two-dimensional' nature  of  the  flow 
can  be  severely  aggravated  as  a  result  of  being  accelerated  by  MHD 
forces . 
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The  comparison  of  theory  with  experiment  in  the  one  case  tested 
reveals  that  the  predicted  impact  pressures  are  greater  on  the  anode 
wall  than  on  the  cathode  wall,  in  accord  with  measurements,  and  the 
average  values  are  comparable.  However,  more  refined  criteria  for 
specifying  the  nonuniformity  in  the  initial  profiles  is  needed  since  there 
is  substantial  disagreement  between  theory  and  experiment  in  the  amount 
of  nonuniformity  in  the  exit  region. 

The  comparison  of  one-  and  two-dimensional  solutions  for  flow 
through  a  linear  Hall  generator  operated  in  a  two-terminal  mode  shows 
that  if  the  upstream  gas  dynamic  state,  seed  rate,  and  terminal  voltage 
are  matched,  the  generated  power  is  larger  in  the  two-dimensional  case 
(16  percent  in  the  test  case  analyzed).  The  two-dimensional  analysis 
gives  the  end  effects  caused  by  (1)  the  fringing  magnetic  field  and 
(2)  the  conduction  of  the  electrical  current  through  the  terminals  into 
the  flowing  gas.  Although  a  small  difference  in  generated  power  was 
experienced  because  of  the  end  effects  in  the  case  considered  (16  per¬ 
cent  as  previously  noted)  each  generator  design  should  be  investigated. 

A  more  important  aspect  of  the  two-dimensional  analysis  is  the  calcu¬ 
lation  of  the  pressure  gradient  transverse  to  the  flow.  This  gradient 
can  have  adverse  effects  on  the  boundary  layer  in  the  channel. 

It  would  not  have  been  possible  to  obtain  any  of  the  two-dimensional 
solutions  presented  in  this  report  without  using  the  direct  method  of 
Ref.  4  to  solve  the  finite  difference  relations.  The  simplicity,  stability, 
and  accuracy  of  the  method  is  such  that  it  is  recommended  in  general 
to  solve  systems  of  second-order  finite  difference  equations  in  prefer¬ 
ence  to  other  methods  that  have  been  proposed  in  the  literature. 
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note:  all  dimensions  in  meters 

OPERATING  GAS  ~  NITROGEN 

SEED  ~  POTASSIUM  (0.50  %  by  Weight) 


a.  Geometry  and  Applied  Fields 
Fig.  1  MHD  Flow  Field  in  Accelerator  A 
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b.  Velocity  and  Temperature;  Uniform  versus  Nonuniform  Initial  Profiles 

Fig.  1  Continued 
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c.  Current  Density;  Uniform  Initial  Profiles 
Fig.  1  Continued 
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d.  Current  Density;  Nonuniform  Initial  Profiles 
Fig.  1  Concluded 
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b.  Velocity  and  Temperature;  Uniform  versus  Nonuniform  Initial  Profiles 

Fig.  2  Continued 
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OPERATING  GAS  —  AIR 
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a.  Geometry  and  Applied  Fields 
Fig.  3  MHD  Flow  Field  in  Accelerator  C 
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b.  Velocity  and  Temperature  without  Finite  Electrode  Segmentation 

Fig.  3  Continued 
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to 

-o 


UPPER  WALL 
(CATHODE) 


<L 


LOWER  WALL 
(ANODE) 


c.  Impact  Pressures  with  Finite  Electrode  Segmentation 
Fig.  3  Continued 
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a.  Geometry  and  Applied  Fields 
Fig.  4  MHD  Flow  Field  in  Hall  Generator 
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w 

o 


ONE- DIMENSIONAL  SOLUTION 


TWO-  DIMENSIONAL  SOLUTION 


p  =  0.39  atm 
T  =  2421  °K 

I 


} 


UPPER  WALL 


p  «  0.45  atm 
T  =  2437  °K 


} 


LOWER  WALL 


p  =  0.204  atm 
T  *  ,2397 


£ 


itm  1 
°K  J 


UPPER  WALL 


Vx  =  10, 000  volts 
I _  -  2316  amps 


S*  =  0.18476  %  by  weight 


p  =  V  237  atm 
T  =  2382  °K 


} 


LOWER  WALL 


NOTE:  GAS  DYNAMIC  STATE  UPSTREAM  OF  THE  MAGNETIC  FIELD  IS  THE  SAME  FOR  BOTH  SOLUTIONS 

b.  One-Dimensional  versus  Two-Dimensional  Solutions 
Fig.  4  Continued 


AEDC-TR-71-1 81 


GROUNDED  SECTION 


POWER  SECTION 


SEGMENTED  ELECTRODES 
Ix a  Generated  Current 


¥  =  O.IIx  Ix  2IX  3IX  4IX  5IX  6  Ix  7  Ix  8IX  8.7  Ix 


c.  Electric  Current  Paths 
Fig.  4  Continued 
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APPENDIX  II 

CONTINUITY  OF  MASS  AND  ELECTRIC  CURRENT 
(SOURCE  FLOW) 


The  usual  forms  in  which  the  mass  and  electric  current  conserva- 

— *  — > 

tion  equations  are  given,  V  •  pV  =  0  and  V  •  j  =  0,  are  applicable  to 
plane  or  three-dimensional  flows.  However,  it  is  incorrect  to  apply 
these  formulas  to  the  flow  model  analyzed  in  this  report  because  the  flow 
is  quasi-two-dimensional  in  the  sense  described  below. 


It  is  intended  that  the  analysis  should  have  the  capability  of  com¬ 
puting  flows  through  channels  with  divergence  of  both  pairs  of  sidewalls 
(the  B-field  as  well  as  E-field  walls)  while  at  the  same  time  expressing 
the  equations  of  motion  as  functions  of  only  two  spatial  coordinates 
(x  and  y).  In  order  to  accomplish  this  it  is  assumed: 

1.  The  x  and  y  components  of  velocity  and  current  density 
are  functions  only  of  x  and  y 

2.  The  ratio  of  the  z  and  x  components  of  velocity  and 
current  density  satisfy 


and 


w  _  z  dW  » 

u  "W  dx  ,''i 

i  ' 

i 

i 

3Z  z  dW  i 


i  i 


This  latter  assumption  states  that  the  slope  of  the  streamline  projection 
in  the  x-z  plane  is  a  linear  function  of  the  z  coordinate. 
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Let  F  =  [Fx,  Fv>  Fz]  represent  either  of  the  two  solenoidal  vectors, 

■4  — >  "4  . 

pV  or  j«  The  combination  of  the  above  assumptions  with  V  -  F(x,  y,  z)  =  0 
yields 


V.F(x,y,z)  =  F  (x, y)  +  ~  F  (x, y)  + 


dx  x 


5 

-  F 

dx  x 


by  y 

b 

+  -  F 

by  y 


a 

dz 

^Fx(x,y) 

z  dW  I 

W  J 

F 

X 

dW 

W 

dx 

1  3 

- WF 

W  bx  x 


3 

+  F 
by  y 


1  ^ 

-  W  WFX 


+  37  WFy> 


The  conclusions  are 


and 


—  Wpu  +  —  Wpv  =  0 
ox  ay 


—  Wj  +  —  Wj  =0 

ax  x  ay  y 


(ii-D 


(II -2) 


Equation  (II- 1)  suggests  the  definition  of  a  (normalized)  mass  stream 
function  which  satisfies 

5*  2W 

- - 7-  pv 

ox  ra 

(II-3) 

b*  2W 

a7'\rpu 


where  the  normalization  has  been  chosen  such  that  ^  =  1  on  the  upper 
wall  and  i>  =  -1  on  the  lower  wall.  Equation  (II - 2)  suggests  the  defini¬ 
tion  of  an  electric  current  stream  function  which  satisfies 


dY(x,y) 

bx 


ll(x,y) 

dy 


(II-4) 
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The  solution  of  the  describing  equations  is  performed  using 
von  Mises  coordinates;  this  requires  to  be  expressed  as  a  function  of 
(x,  ip)  as  follows 

dY(x,i|r)  _  [dx(s,n)|  w . 

^  "Us  J 

<3  ^  J  X 

where  (s,  n)  are  orthonormal  coordinates  measured  tangential  and 
normal  to  fluid  particle  paths. 
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APPENDIX  III 

TRANSFORMATION  OF  EQUATIONS  OF  MOTION  , 


The  gas  dynamic  equations  are  solved  after  rearranging  the  steady- 
state  form  of  the  conservation  equations  for 

Mass  (Source  Flow) 


Momentum 


Energy 


7-  WpV  =  0  « 


„  DV  -3 

p  —  +  Vp  =  JXB 

Dt 


(III-l) 


(III -2) 


(HI -3) 


into  a  system  of  differential  equations  where  the  dependent  variables 
are  q(x,  <p),  p(x,  ip),  T(x,  ip)  in  the  accelerator  analysis,  and  q(x,  <P), 
p(x,  ip),  H(x,  ip)  in  the  generator  analysis.  These  variables  were 
selected  because  they  are  of  special  importance,  and  because  the 
thermodynamic  functions  are  expressed  as  functions  of  p  and  T  in  the 
accelerator  analysis,  and  p  and  h  in  the  generator  analysis. 


The  equation  for  q  is  obtained  by  combining  the  components  of  the 
momentum  Eq.  (Ill -2)  to  yield  (in  von  Mises  coordinates) 

d  q2  dp  , 

pu  —  —  +  u  =  [VjB]  (III-4) 

dx  2  dx 

The  assumption  that  B  has  a  component  in  the  z -direction  only  implies 


where 


[vjsi  =  (|2)'1  uJnB 
dx  .  dy.  . 

Jn  ds  Jy  "  ds  Jx 


(HI-5) 


dx 

ds 


dy 

ds 


dy 

35 


{■  *  <!>' } 


1/2 
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and 


^yCxjifr)  is  given  by  Eq.  (4). 

dx 


The  result  is 

^q(x,f) 

dx 


(pq)-1 


[ 


os 


bp] 
bx  i 


in  which  form  the  solution  is  obtained  using  a  Runge-Kutta  method,  and 
the  pressure  gradient  that  appears  on  the  right-hand  side  is  evaluated 
below. 


The  pressure  is  computed  on  two  levels  of  approximation.  The  first 
level  of  approximation  consists  of  solving  a  system  of  streamline  differ¬ 
ential  equations  for  pressure.  This  is  followed  by  a  second  level  of 
approximation  in  which  the  pressure  is  recomputed  by  integrating  the 
transverse  momentum  equation  in  the  y  direction.  These  approxima¬ 
tions  are  obtained  from  the  following  sequence  of  calculations:  con¬ 
servation  of  mass  is  expressed  in  integral  form;  this  integral  relation¬ 
ship  is  differentiated  in  the  streamline  direction;  by  invoking  an  equation 
of  state,  and  momentum  and  energy  conservation,  a  system  of  integral 
equations  for  streamline  pressure  gradient  is  obtained;  these  integral 
equations  are  then  formally  inverted  to  obtain  the  required  formula  for 
pressure  prediction  to  be  performed  in  an  explicit  manner. 


The  detailed  derivation  of  the  formula  for  ^  p^(x,  p),  where  p^  is 
the  pressure  on  the  ith  streamline,  begins  by  differentiating 


The  next  step  is  to  substitute  into  the  previous  equation  either 
Eqs.  (Ill - 6)  to  (III- 10)  in  the  case  of  the  accelerator  analysis. 
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P  =  p(P,T) 

( HI  -  6) 

Sp  _  Sp  dp  S£  ST 

Sx  Sp  Sx  St  Sx 

(III  -  7) 

St  _  (Sh,-1  f(I  .  Sh)  S£  +  5-7  -  IvjB]  ]  (m_8) 

Sx  St  L  p  ^p  pu  J 

(Eq.  ( III  -  8 )  follows  from  Eqs.  (8)  and  (III -4)) 

..2  2  TSp  Sp  Sh.  Sh  -li 

M  [a;  +  5  (1  -  p  a?' (p  a?  J 

( III  -  9) 

Fi  -  (— )~*  1e 

St  St 

(III-10) 

or  Eqs.  (Ill -11)  to  (III - 15)  in  the  case  of  the  generator  analysis 

p  -  p(P,h) 

(III- 11) 

Sp  Sp  Sp  +  Sp  Sh 

.  Sx  Sp  Sx  Sh  Sx 

(m-i2) 

Sh  "j  S  q2 

Sx  pu  Sx  2 

(Eq.  (Ill -13)  follows  from  Eq.  (8)) 

(III -13) 

Ma  =  SL-  fi.pl.  (RZD  -  —  (EZT) 
RZT  |_  Sp  Sh 

(III- 14) 

F 

Fi  =  - 

Sh 

In  addition  both  analyses  use  Eq.  (16),  which  implies 

!*•  J  =  E- j  -  [v j B ] 

(III-15) 

Jo'  +  3 


n 
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and  the  identity  which  makes  it  possible  to  formally  invert  the  pressure 
integral  equations 


3 

—  P. 
dx  J 


P. 

i 


to  achieve  the  result 


2  dA 

-T“  +  F2-Fa 

in  dx 


—  P,  = 
dx 


pvi'1  « 


(M  -  1)  d* 


where 


-1  ds 
1 

Fs  -  f  (pq)  ~ 1  ~  (■—)  d\jr 


dx  ds 


1 'vr'{ 


(MS  -  1)  (Pi 
dx  J 


.dx 


- 1 


+  (^)  (jnB  +  P'lq  E*-j  Fi))-d*3 


>}' 


The  notation  used  for  the  dummy  variable  of  integration  in  the  second 
integral  of  the  numerator  means  that  the  quadrature  is  to  be  performed 
with  the  streamlines  labeled  with  j  subscripts.  This  is  done  in  order  to 
define  the  meaning  of  pj  -  pj  which  occurs  in  this  same  integral;  thus 

Pj  -  pj  means  the  difference  between  pressures  on  the  jth  and  ith  stream¬ 
lines.  Formally,  this  pressure  difference  can  be  evaluated  from  the 
transverse  momentum  equation  where 

dv  di|j  dp 

^  dx  dy  df  x 


is  integrated  in  the  transverse  direction  and  then  differentiated  in  the 
axial  direction  to  give 
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3  v 

and  -5—  is  evaluated  from 
9x 


-1/2 


q 


9y 

and  the  streamline  slope  is  obtained  from  Eq.  (4).  For  channels 

whose  conductor  walls  are  only  slightly  diverged  such  as  that  of  Ref.  11, 
the  last  term  on  the  right-hand  side  of  Eq.  (Ill- 16)  =  0. 


The  second  approximation  for  pressure  alluded  to  above  is  given  by 
the  integral  of  the  radial  momentum  equation 


The  difference  of  current  fluxes  on  the  ith  and  jth  streamlines,  -  Vj, 
is  obtained  in  the  case  of  Faraday  channels  by  the  theory  described  in 
Appendix  V,  and  from 


for  channels  in  which  the  describing  electric  equation  is  that  of  an  elec¬ 
tric  potential. 


40 


A'ED  C -T  R  -7 1  -  f  8 1 


APPENDIX  IV 

THERMODYNAMIC  FUNCTIONS 


Expressions  for  thermodynamic  functions  of  air  and  nitrogen, 
respectively,  are  derived  in  closed  form  in  such  a  manner  as  to  ap¬ 
proximate  their  tabulated  values  given  in  Refs.  6  and  7.  The  starting 
point  in  the  derivation  is  to  develop  formulas  relating  the  thermodynamic 
variables  p,  T,  p,  Z,  h,  and  S  which-  will  ( 1)  approximate  the  tabulated 
values  and  (2)  satisfy  the  thermodynamic  relationships 

p  =  RZpT 

and 

T  dS  =  dh  -  — 

P 


It  facilitates  the  development  of  the  results  to  combine  the  previous 
two  equations  with-  the  postulated  relationships 


|  =  f  i(T,  Z) 
p'  =  f2(T,Z) 

and  the'  identities 

dS  dS 

dS  =  —  dT  +  — —  dZ 
dT  dz 

and 

d  ds  d  ds 
dz  =  d? 


(IV-1) 


which  yields 

=  J_  A  il 

dT  T2  dz  R 


(IV-2) 


as  a  differential  relationship  which  must  be  satisfied  by  the  assumed 
empirical  relationship.  At  this  point  the  remaining  details  of  the 
derivation  become  peculiar  to  the  particular  gas  being  considered. 


AIR  THERMODYNAMIC  FUNCTIONS 

Let  the  static  enthalpy  have  the  functional  form 

£  =  f (T)  +  (Z  -  1)  g(T) 

K 


41 


AEDC-TR-71-1 81 


When  this  is  substituted  into  Eq.  (IV- 2)  the  integrated  result  is 


P 


F(Z>  e 


dT 


The  arbitrary  functions  appearing  in  the  two  previous  equations  are 
evaluated  by  matching  the  exact  values  in  the  higher  range  of  tempera¬ 
tures  representative  of  MHD  operation.  However,  the  results  are  in¬ 
sufficiently  accurate  at  lower  temperature  values.  In  order  to  extend 
the  range  of  validity  a  separate  empirical  fit  is  performed  to  yield 
comparable  low  temperature  accuracy  while  at  the  same  time  render¬ 
ing  all  the  thermodynamic  functions  continuous  ( as  functions  of  both  p 
and  T)  at  a  matching  temperature.  For  the  sake  of  simplicity  the  low 
temperature  formula  for  enthalpy  (that  is  used  when  T  <  Tmatch)  is  not 
required  to  satisfy  the  entropy  condition  of  integrability  (Eq.  (IV- 1)). 
Thus,  the  solution  of  the  equations  of  motion  for  an  adiabatic,  reversible, 
flow  when  obtained  using  these  functions  is  exactly  isentropic  when 
T  >  Tmatch,  but  only  approximately  so  (typically  in  error  by  1  percent) 
when  T  <  Tmatch.  The  results  are: 


T 

match 


4.05  x  104 


56,500 

(1.2  -  Z)  e - T 

(Z  -  l)2 


T  >  T 

match 

^  =  4.96T  +  (Z  -  1) (56, 500  +  T)  -  2000 
|  =  3.96  In  T  +  (Z  -  1)(1  +  — >»5-00-)  -  In  p 
-  0.2  In  (1.2  -  Z)  -  0.919146 


T  < 


T 

match 


j  -  (ai  +  a2T)T  +  (Z  -  1)(56,500  +  T)  +  aa 


|  =  (ai  -  1)  In  T  +  2  a  2T  +  (Z  -  1)  (1  +  — 
-  In  —  -  0.2  In  (1.2  -  Z)  +  a4 

P» 

00 
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Solve  for  Tmatch  from 


a i  +  2a2  T 


match 


4.96 


a i  =  3.322661 

as  =  3. 24611xl0“4 

aa  =  64.6859 
a4  -  9.9744 

The  values  of  and  a.2  are  determined  from  the  observation  that 

h. 

at  sufficiently  low  temperatures  versus  T  is  virtually  a  linear  func- 
'h 

tion;  i.  e. ,  =  aj  +  a2T.  The  value  of  Tmatch  is  determined  by  re¬ 

quiring  continuity  in  Cp;  a.$  and  are  given  by  fulfilling  continuity  of 
h/R  and  S/R,  respectively,  at  the  matching  temperatures. 


NITROGEN  THERMODYNAMIC  FUNCTIONS 


A  single  set  of  expressions  is  obtained  for  the  entire  range  of 
interest  by  postulating  the  expressions  for  h/R  and  Z  given  below  where 
the  constants  are  evaluated  by  optimizing  the  accuracy  of  Cp. 


h 

R 


Ci 


(e 


6v 

7T 


-  1)T 


I 


+  C2T  +  Z  T  +  Ca  (Z  -  1) 


Z-l 


Cg 


LC4  T  J 


c7 


=  (1  + 

Cl)  In  T  +  2C2T 

-  in  £- 

P. 

+  (Z  ■ 

jl 

-  mi  +  -  +  _> 

e 

e 

6 

+ 

v 

—  coth 

°v 

v  v 

ln(sinh  — )  +  Ca 

2T  ■ 

2T 

2T 

Cl  = 

2.5 

C0  = 

.496404 

c2  = 

7. 57143xl0-8 

c7  = 

55830. 65 

c3  = 

c7/c8 

C8  = 

2.38 

c4  = 

3. 66099xl0“3 

1! 

> 

ct> 

3340 

Cb  = 

?76. 6027 
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APPENDIX  V 

ELECTRIC  CURRENT  DENSITY  IN  FARADAY  CHANNELS 


Since  the  boundary  conditions  for  the  electrical  portion  of  the  prob¬ 
lem  of  flow  through  Faraday  channels  is  posed  in  terms  of  current  (on 
the  insulator  portion  of  the  walls  jn  =  0,  along  the  electrode  section 
jn(x)  is  specified,  and  at  the  upstream  and  downstream  ends  of  the 

region  of  computation  j  =  0),  it  is  natural  to  solve  a  differential  equa¬ 
tion  for  a  current  stream  function  (rather  than  an  electric  potential) 
since  this  results  in  the  boundary  conditions  being  of  the  Dirichlet  type. 
The  derivation  of  this  equation  is  performed  by  substituting  E  as  a 
function  of  j,  from  Ohm's  law,  into  the  steady-state  form  of  Faraday's 
law,  and  then  substituting  a  current  stream  function  for  j  from  con¬ 
servation  of  current.  The  analytical  expressions  of  these  statements 
are 


V  x  E  =0 

aEs  "  4  +  ,JJT  ^ 

a(En-qB)  =  -  cot  +  jn 

Wj  «  V  x  [0,  0, 


(V-l) 


Since  the  curl  operation  is  invariant  under 
coordinate  transformation  it  follows  that 


df (s , n) 


an  orthonormal  curvilinear 


and 


5Y(s,n) 

ds 


When  a  further  coordinate  transformation  to  the  von  Mises  variables, 
with  which  the  problem  is  solved,  is  introduced  the  result  is 


and 


,  dy  5  df  3 .  . 

WJS  =  (_  7"  7"  +  T*  T')  fCx’,lr) 

3  os  ox  dn  oilr 


nx’*’ 


( V-2) 
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The  combination  of  Eqs.  (V-l)  and  (V-2)  results  in  the  elliptic  differen¬ 
tial  equation 


where 


¥  + 
xx 


+  A 


3  Y  + 
X 


=  ae 


+ 


Ao  = 


(5i> 

dn 


A3 


=  -  uut  ( — 

L  as 


■  a2y  ^ Wa  w 

)  — —  +  -  +  Wa 

dx2  Wa 


diif  a  c 

dn  S\|r  L 


bx  by 

(JOT  -5—  +  -ri¬ 
ds  os 

Wo 


}] 


A4  -  Wa 


[  bx  by  b 

(a,T  55  '  ^  55 


dt-  dir 

,  cfn  ^  d*  °  ,  dn  v 

( - )  +  —  —  ( - ) 

Wa  dn  diji  Wa 


dx  aujr  dO 

3s  bx  dn 


As  "  w°  If  s qB 


As  -2^®* 
ds  dn 

The  term  containing  the  mixed  derivative  is  included  on  the  right-hand 
side  where  it  can  be  treated  as  a  known  quantity  in  the  method  of 
successive  approximations  that  is  used  to  simultaneously  solve  the 
coupled  gas  dynamic -electrical  problem.  The  motivation  for  proceed¬ 
ing  in  this  manner  is  that  the  inclusion  of  mixed  derivative  terms  on  the 
left-hand  side  adds  complexity  to  the  finite  difference  method  of  solu¬ 
tion. 


45 


AEDC-TR-71-1 81 


APPENDIX  VI 

ELECTRIC  POTENTIAL  IN  HALL  CHANNELS 


The  boundary  conditions  for  the  electrical  portion  of  the  problem  of 
flow  through  the  Hall  channel  considered  in  this  report  (two-terminal 
operation  with  equipotential  surfaces  on  the  ends: 

=  j  (x)  . 

upper  n  lower 
wall  wall 


V  =0  and  j  (x) 
y  n 


in  the  electrode  region)  are  such  that  it  is  natural  to  solve  a  differen¬ 
tial  equation  for  an  electric  potential  (rather  than  for  a  current  stream 
function)  since  this  results  in  the  boundary  conditions  on  the  ends  being 

of  the  Dirichlet  type.  The  derivation  of  this  equation  is  performed  by 
— ►  — ^ 

substituting  j  as  a  function  of  E,  from  Ohm's  law,  into  the  equation 
describing  conservation  of  current,  and  then  substituting  an  electric 
potential  for  E  as  a  consequence  of  Faraday's  law.  The  analytical  ex¬ 
pressions  of  these  statements  are 

V.Wj  = 

■^s 


E  = 

Since  the  gradient  operation  is  invariant  under  an  orthonormal  curvi¬ 
linear  coordinate  transformation,  it  follows  that 

dcp(s,n) 

Es  "  ‘  * 

5cp(s,n) 

E  =  -  — 

dn 


- - -  k  ~  u>t(e_  -  qB)l 

1  +  (u>t)2  Is  n  J 

[E  -  qB  +  hit  E 
n  sj 


1  +  (out) 
-  vcp 


(VI-1) 


When  a  further  coordinate  transformation  to  the  von  Mises  variables  is 
introduced,  the  result  is 


E 

s 


3x  3 

(-  r-  T”)  V(*,0) 
os  ox 


(VI-2) 


E 


n 


ds  dx 
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The  combination  of  Eqs.  (VI- 1)  and  (VI-2)  results  in  the  elliptic 
differential  equation 


where 

A  i  = 


+  AB«PH  +  AsCPx  +  ^  =  As  +  As  cp^ 
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Wa 


dil/  2 
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APPENDIX  VII 

DIRECT  METHOD  OF  SOLUTION  OF  FINITE  DIFFERENCE  EQUATIONS 


The  partial  differential  equations  for  electric  current  stream  func¬ 
tion  and  potential  that  are  derived  in  Appendixes  V  and  VI  are  written  in 
implicit  finite  difference  forms  using  central  differences  as  approxima¬ 
tions  to  the  derivatives.  The  resulting  matrix  equation  is  block  tridiag¬ 
onal;  i.  e. ,  the  partitioned  form  of  the  coefficient  matrix  is  tridiagonal. 
The  methods  that  have  been  proposed  in  the  literature  for  solving  such 
an  equation  (i.  e.  ,  a  system  of  linear  algebraic  relationships)  fall  into 
the  general  categories  of  iterative,  direct,  and  block  iterative  which  is 
a  hybrid  of  the  first  two  methods.  Several  iterative  schemes  have  been 
applied  by  various  authors  to  the  solution  of  implicit  finite  difference 
equations  when  their  numerical  properties  were  sufficiently  simple.  An 
extensive  test  of  block  iterative  techniques  in  the  course  of  the  analysis 
that  forms  the  subject  of  this  report  was  made  and  it  was  found  that 
iterative  methods  are  incapable  of  solving  the  electrical  boundary  value 
problems  that  are  described  herein.  The  problems  were  ultimately 
solved  by  the  use  of  a  direct  method  on  a  numerical  gridwork  whose 
size  is  limited  only  by  the  available  memory  of  the  digital  computer 
being  used.  It  has  been  shown  by  substitution  that  the  values  computed 
for  the  unknown  variable  invariably  satisfy  the  given  matrix  equation. 

That  the  success  of  such  an  approach  was  not  a  priori  obvious  is  ex¬ 
plained  by  the  well-known  fact  that  the  finite  word  length  of  a  computer 
results  in  an  accumulation  of  truncation  error  that,  in  general,  pre¬ 
cludes  the  successful  employment  of  a  direct  method.  However,  it  was 
found  that  one  particular  variant  of  the  direct  algorithm  proposed  by 
Schechter  in  Ref.  4  has  the  capability  of  negating  truncation  error  accumu¬ 
lation.  The  method  proceeds  by  accomplishing  a  forward  elimination 
with  a  conventional  factorization  of  the  coefficient  matrix  into  a  product 
of  lower  and  upper  diagonal  matrices.  The  distinguishing  character¬ 
istic  of  the  Schechter  method  is  that  instead  of  performing  the  entire 
backsweep  by  obvious  recurrence  relations  (w’hich  would  invariably  lead 
to  numerical  instabilities),  the  backsweep  operation  is  periodically 
interspersed  with  formulas  for  the  unknowns  that  involve  the  use  of 
intermediate  calculations  performed  and  stored  in  the  computer  memory 
during  the  forward  elimination  process.  Thus  by  limiting  the  number  of 
consecutive  applications  of  straightforward  recurrence  relationships,  it 
has  been  found  possible  in  all  cases  to  perform  the  calculations  within 
numerical  stability  boundaries.  The  efficacy  of  the  method  is  attribut¬ 
able  to  this  property  and  also  to  the  simplicity  of  the  algorithm,  which 
is  presented  below  in  the  notation  of  Ref.  4.  In  order  to  display  this 
simplicity  the  entire  matrix  algebra  necessary  to  effect  the  solution 
{given  by  Eqs.  { VII - 1 )  through  (VII-5))  is  presented  first,  then  a  defini¬ 
tion  of  terms  is  given. 
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FORWARD  ELIMINATION 


Ai  =  Mi 


A  s  M 
n  n 


W"H 

{yn}*{gn} 


-D  A  -1  E  ,,  1  <  n  :£  q,  (VII- 1) 

n  n-1  n-1  ’  k 


-  D  A  _"1  { y  1  <  n  £ 

n  n-1  1/n-lJ  ’ 


(VII -2) 


BACKSWEEP 


{\kY  V1  {%} 


fn-J-  -  Mn  {%}  *  E„  K+l}}  <Vn- 


(VII -3) 

4) 


{vn}“  V1  K}  -  V1  Er  {vn+l} 


(VII -5) 


The  braces  around  the  terms  indicate  vector  matrices  of  length  p,  and 
the  other  matrices  are  rectangular. 


The  symbolism  used  in  Eqs.  (VII- 1)  through  (VII-5)  (which  is  the 
same  as  that  used  in  Ref.  4)  is  related  to  the  physical  problem  as 
follows.  Denote  the  first  upstream  column  of  gridpoints  by  n  =  1  and 
the  last  downstream  column  by  n  =  where  k  denotes  the  number  of 
changes  in  the  number  of  unknowns  to  be  solved  for  on  the  individual 
columns.  Let  pr  denote  the  number  of  rows  on  column  qr._ 


_L 

r 

rows 

L 

1  i 

n=l  n=qi 


r 

pr  rows 

_ L_ 


i  i 

n=q±  n=qfc 


Physical  Array 

As  is  well  known  the  finite  difference  relations  for  a  linear  second - 
order  partial  differential  equation  without  mixed  derivatives  can  be 
written  in  partitioned  matrix  form  as 
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where  Dn,  Mn,  En  are  matrices  with  the  same  number  of  rows, 

En,  Mn+1,  Dn+2  have  the  same  number  of  columns,  and  the  Mn  are 
square  (pr  x  pr).  The  vectors  |vn  j  are  the  unknowns;  in  the  context  of 

this  report  these  are  successive  approximations  to  the  current  stream 
function  (Faraday  channels)  or  electrical  potential  (Hall  channels). 


The  derivation  of  the  solution  algorithm  is  performed  as  follows. 
Denote  the  coefficient  matrix  by  Q;  Eq.  (VII -6)  can  then  be  written  as 

Q  {v}  =  jgj  .  Equation  ( VII - 7 )  results  from  the  factorization  of  the 

coefficient  matrix  Q  into  a  product  of  lower  and  upper  diagonal  matrices 
Q  =  LU. 
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(VII -7) 


Equation  (VII-2)  is  obtained  by  solving  for  |y|  from  |y  J  =  U  |vj  and 
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g 

Equation  ( VII - 3 )  is  obtained  from  the  last  row  of 
Equation  (VII-4)  is  obtained  by  rewriting  the  nth  row  of 
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Dn  {v„-l}  +  M„  {%}  ♦  E„  {%+1}  -  {*„} 

Equation  (VII-5)  is  obtained  by  solving  for  |vnj  from  the  nth  row  of 

u  f  v  }  =  fy  }.  which  is  An  |vn  }  +  En  {vn+1  }  =  |ynJ.  The  reader  is 

referred  to  Ref.  4  for  additional  details  of  the  matrix  algebra  that 
are  required. 
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APPENDIX  VIII 


TABLE  I 

INITIAL  CONDITIONS— UNI  FORM  AND  NONUNIFORM  INLET  PROFILES, 

ACCELERATOR  A 


The  initial  conditions  listed  below  are  values  at  the  entrance  to  the  electrode  section 
that  would  result  from  an  isentropic  expansion  through  the  channel. 


Ul 

to 


UNIFORM  INITIAL  PROFILES 


p/p  =0.5 

09 

M  =  1.6 


Assumptions  ^ 


H  =  5.4168  (10®) 
m2/sec2 


NONUNIFORM  INITIAL  PROFILES 


p/p„ 

00 


=  0.5 


Assumptions  < 


M  =  1.6 

M  is  uniform  in 

transverse  direction 

T  is  a  quadratic 
function  of  the 
transverse  coordinate 


"  1*5  X  Twall 
T^  -  3600#K 


The  assumptions 
m2/sec  ,  and 


imply  H 
a  veloc?¥f 


=  5.4168  (10®) 
profile. 


The  assumptions  imply  values  of 
velocity  and  temperature. 
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TABLE  II 

INITIAL  CONDITIONS— UNIFORM  AND  NONUNIFORM  INLET  PROFILES, 

ACCELERATOR  B 


The  initial  conditions  listed  below  are  values  at  the  entrance  to  the  electrode 
section  that  would  result  from  an  isentropic  expansion  through  the  channel. 


UNIFORM  INITIAL  PROFILES 


NONUNIFORM  INITIAL  PROFILES 


f  P/P.  =  °-5 


p/p  -  0.5 

GD 


cn 

oo 


Assumptions  i 


u  =  3000  m/sec 


T  =  3400°K 


The  assumptions  imply 


Assumptions  ^ 


p  is  uniform  in 

transverse  direction 

H  is  a  quadratic  function 
of  the  transverse 
coordinate 


m  =  0. 7158  kg/sec 

H  «  10.262  (10® ) 
m2/sec2 


Hg  1.5  x  Hwall 
m  =  0.7158  kg/sec 
Hayg  -  10.262  (106)  m2/sec2 

b 

The  assumptions  imply  profiles  of  velocity 
and  temperature. 
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TABLE  III 

INITIAL  CONDITIONS— NONUNIFORM  INLET  PROFILES, 

ACCELERATOR  C 

The  initial  conditions  listed  below  are  values  at  the  entrance  to  the  electrode  section 
that  would  result  from  an  isentropic  expansion  through  the  channel. 

The  following  one-dimensional  values  are  from  Ref.  11,  Run  1323: 

p/pM  =  4.7 


T,  =  3100°K 

1-dim 

u  =  2960  m/sec 

1-dim 

The  criterion  for  choosing  the  following  relationships  is  explained  in  Ref.  5. 
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